Computational evaluation of zirconocene catalysts for ε-caprolactone cationic ring-opening polymerization

This quantum chemical study presents the ligand effect and a structure–property relationship in the cationic ring-opening polymerization (CROP) of ε-caprolactone using zirconocene catalysts. We first examined the effects of catalyst structure on the initiation and chain propagation steps of the CROP process. A total of 54 catalyst structures were investigated to understand the influence of the ligand structure on the stability of the catalyst–monomer complex and polymerization activity. The properties of the catalysts were analyzed in terms of ancillary ligands, ligand substituents, and bridging units. Calculations showed that the polymerization follows a proposed cationic mechanism, with ring opening occurring via alkyl-bond cleavage. A correlation between complex stability and activation energy was also observed, with ligand substituents dominating in both steps. While the ancillary ligands directly affect the HOMO energy level, the bridges are mainly responsible for the catalyst geometries, resulting in reduced complex stability and higher activation energy for the propagation step. This study contributes to a better understanding of the structural characteristics of zirconocene catalysts, which offers guidance for improving CROP activities in lactone polymerization.

www.nature.com/scientificreports/polymerization of CL in toluene at room temperature or 60 °C was obtained with afforded monodisperse PCLs (M w /M n = 1.06-1.3).The best catalytic system was found for the methyl substituted ligand, (Me 5 Cp)CpZrMe 2 , as it polymerized CL with good control (M w /M n < 1.13, M n,calcd ≈ M n,GPC ), and the molecular weight of the resulting PCLs increased with monomer conversion.They also suggested the higher catalyst efficiency as (Me 5 Cp) 2 ZrMe 2 > (Me 5 Cp)CpZrMe 2 > Cp 2 ZrMe 2 .Ten years later, a research group of Pitsikalis et al. 24 reported the CROP of CL and δ-valerolactone using zirconocene catalytic systems with three different borate cocatalysts, B(C 6 F 5 ) 3, [B(C 6 F 5 ) 4 ] -[Ph 3 C] + and [B(C 6 F 5 ) 4 ] -[Me 2 NHPh] + .In this work, three zirconocenes were introduced in the reaction, namely the bis(cyclopentadiene) ligand (Cp 2 ZrMe 2 ), bis(indene) ligand (Ind 2 ZrMe 2 ), and tertiary butyl substituted ligand (( t BuCp) 2 ZrMe 2 ) and the high molecular weights PCL and narrow molecular weight distributions (M n = 10,000-38,000 g/mol and M w /M n < 1.1) were obtained in a well-controlled manner.The catalyst performance was also reported as Cp 2 ZrMe 2 < Ind 2 ZrMe 2 , and ( t BuCp) 2 ZrMe 2 .The same research group also applied their catalytic systems for copolymerizations of other monomers such as methacrylate 24 , vinyl ether 32 and 2-oxazoline 23,31 .The first mixed cyclopentadienyl-indenyl zirconocene, CpIndZrMe 2 , was used as a sole catalyst for the polymerizations of CL without adding any cocatalyst 33 .The resulting polymers yielded high quantitative outcome in terms of medium molecular weight polymers and moderate to broad polydispersities (M n = 25,000-128,000 g/mol and M w /M n = 1.11-1.98).These experimental studies open the way to enhance catalyst activities, polymer yields and molecular weights of the zirconocene polymerization by ligand modifications.
A number of density functional theory (DFT) studies on the metallocene-initiated CROP of cyclic esters and related monomers have also been reported (see references therein 34,35 ), which provided valuable aspects of DFT modeling and visualization of cationic mechanisms and molecular structures in complement to experimental observations.For example, the electronic and geometric properties of dimethylzirconocene catalysts were greatly influenced by the presence of the silicon-bridge and ancillary ligands 19 .The effects of ligand structure, solvent, and metal on the dimerization of dinuclear zirconocene cations have also been investigated with DFT 36 .We also studied the influence of ligand structure on the initiation and propagation steps of the trimethylene carbonate polymerization using the naked model approach 17 .The CROP process follows an active chain end (ACE) mechanism where propagating species are tertiary oxonium ions located at the chain end (Fig. 2) 37,38 .In brief, a precatalyst Cp 2 ZrMe 2 is activated by a borate cocatalyst (e.  to the initiation (or complexation, Cation → Complex) and propagation steps (Complex → R1 → P1), depicted in Fig. 2.
Here, we employed the DFT modeling approach to study the effects of ligand structure through a series of zirconocene complexes (varying in steric hindrance and electronic properties) on the initiation and propagation.Their interactions with the boron cocatalysts on the rate-limiting step are analyzed in terms of the relative complex stability and the activation energies of the process.The structure-property relationship was then performed for the 54 zirconocene structures and the results were interpreted in terms of ancillary ligand, ligand substituents and bridge units.This study provides a better understanding of the structural characteristics of zirconocene catalysts in CL polymerizations.

Computational methods
All geometries were optimized without restraints in gas-phase using density functional theory (DFT) 39,40 at the B3LYP level.This B3LYP is chosen because of its success in DFT modeling studies on similar processes 7,13,17,35,41 and its lower activation energy compared to other DFT functionals (Figure S1 of Supporting Information, SI).Zirconium atoms were treated with the effective core potential double-ζ basis set (LANL2DZ) 42 whereas all non-metal atoms (C, H, O, B, and F) with a double-ζ basis set, 6-31G(d).The DFT/mixed basis set method has been largely applied for transition metal complexes 12,13,18,19,43 .Harmonic frequency calculations were performed to confirm the transition states and local minima obtained and to determine the Gibbs free energy (G) at standard temperature and pressure (298.15K and 101.325 kPa).Accurate free energies and molecular electronic properties (Mulliken population analysis, HOMO-LUMO energies, and dipole moment) were obtained at the CPCM(toluene)-B3LYP-D3(BJ)/6-311++G(d,p) (LANL2DZ) level using single-point calculations on the optimized structures, incorporating implicit solvation and dispersion correction.All calculations were carried out using the Gaussian 09/16 software 29, 44 .

Results and discussion
In this study, we evaluated the effects of zirconocene structures on the process by means of DFT calculations.Initially, we considered five catalysts that have been previously reported for their CROP activities 24,25 and dimeric properties 36 , namely Cp 2 ZrMe 2 (C1), (Me 5 Cp)CpZrMe 2 (C3), (Me 5 Cp) 2 ZrMe 2 (H1), (tBuCp) 2 ZrMe 2 (C4) and Ind 2 ZrMe 2 (H16) (Fig. 3).These catalysts were selected for comparison purpose, and their properties vary in steric hindrance (C1, C3, H1, and C4) and electronic properties (C1 and H16).Their interactions with the two borate cocatalysts, [MeB(C 6 F 5 ) 3 ] -and [B(C 6 F 5 ) 4 ] -, in the initiation and propagation were analyzed for the activity of the process.Detailed information on how to obtain the structures of zirconocene/borate systems can be found in SI.

Catalyst structure influence on the process in the presence of cocatalyst
To evaluate the effects of zirconocene structures on the initiation and propagation using zirconocene/ borate catalytic systems, five catalysts (C1, C3, H1, C4, and H16) that have been previously used in lactone polymerization 24,25 were considered (see Fig. 3).Calculations of the initiation and propagation mechanisms were carried out using a similar procedure as previously reported in the literature 17 .Values of complexation energies for the Cation → Complex step (ΔG com ) and activation Gibbs energies for the first and second propagations (ΔG ǂ Ea1 and ΔG ǂ Ea2 ) using the five zirconocenes of three different catalytic systems (C + , 1 and 2) are also included in Table 1.The heatmap corresponding to the energy cross-correlation between the catalyst pair for the Figure 3. Five zirconocenes included for the study of catalyst structure effect.Nomenclatures used for each catalyst (C1, C3, H1, C4 and H16) were assigned according to our recent DFT study on the zirconocene dimer properties 36 .The order of the dimer stability of these catalysts is C1 > C3 > H16 > > H1 > C4.
Vol:.( 1234567890 The ion-pair interactions between the different catalyst cations and the cocatalyst anions, [MeB(C 6 F 5 ) 3 ] − and [B(C 6 F 5 ) 4 ] − , during the reaction were analyzed in terms of relative complex stabilities, activation free energies, and geometries, in comparison with the available experimental data.
For the chain propagation, the activation free energies generally follow a trend C + > 1 ~ 2 (Table 1 and Fig. 4b): the calculated barriers, ΔG ǂ Ea1 , for C1, C3, H1, C4, and H16 ranged from 74.0-92.5 kJ/mol for C + to 71.6-88.9 and 70.6-83.2kJ/mol for 1 and 2, respectively.This effect is more profound for the second propagation (ΔG ǂ Ea2 : and positive ΔΔG, Table 1; Figure S4).The bar charts in Fig. 4c further show that the relative barrier heights, ΔG ǂ Ea1 , follow the order C1 < C3 < H1 for Cp 2 ligand and C1 < H16 for Ind 2 ligand, which account for the activity difference (e.g., M n , M w /M n ; Table S1) reported by Hayakawa et al. 25 and Pitsikalis et al. 24 .These results emphasize the influence of steric hindrance and bulkiness of the catalyst structures, which can be used as an important factor for ligand modifications.

Geometric impact
The optimized transition-state structures (1-TS1 and 2-TS1) in the first propagation step of 2) systems (X = C1, C3, H1, C4, and H16) are depicted in Figs. 5 and 6, respectively.Important geometric parameters for the first and second transition-state structures (TS1 and TS2) for the three catalytic systems (C + , 1, and 2) are included in Table 2.As expected, the ligand structure affects mainly the geometries around the catalytic center of the propagating chain end, with small changes in the C α -O α distance among the C + , 1, and 2 systems.The presence of the counteranion results in the elongation of the C α -O α bond by 0.02 Å, compared to that in the C + system, and this bond becomes longer at TS1, with increasing of the methyl group on the Cp ligand (e.g., 2.049 Å at 1-TS1 C1 , 2.078 Å at 1-TS1 C3 , and 2.082 Å at 1-TS1 H1 ).Similar trend was also found in the system 2 (Fig. 6).This is in agreement with the observed difference in the catalyst activity 25 .
The steric effect on the ligand structure is greater in C4, which has steric crowding on its Cp ligand.This limits a cocatalyst access near the C α atom of the propagating chain end, as evidenced by the longest C α -B distance of 7.2 Å (see 1-TS1 C4 in Fig. 5), compared to that in other systems.

Table 1.
Values of complexation energies (ΔG com , Cation → Complex) and the activation energies in the first and second propagations (ΔG ǂ Ea1 and ΔG ǂ Ea2 ), as calculated for the initiation and the first and second chain propagation steps for zirconocene cation, X (C + ), X/MeB(C 6 F 5 ) 3 ˉ (1) and X/B(C 6 F 5 ) 4 -( 2) systems (X = C1, C3, H1, C4 and H16).All energies are in units of kJ/mol, while those in parenthesis are the energy for each system relative to the energy of its counterpart C1.

com/scientificreports/
There is a clear trend between the Zr-B distances and the steric degree on the Cp′ ligand, with Zr-B being significantly increased in the TS2 for bulkier Cp′ ligands, i.e., H1, C4, and H16.For example, the catalyst H1, bearing a highly steric congestion on both Cp′ rings, yields the highest Zr-B distances of ~ 8.5 to 9.5 Å for both 1 and 2 (see Zr-B in Table 2).At 1-TS1, Zr-B and B-C α distances decrease in the order C1 > H1 > H16 > C3.Overall, the results show that the size and shape of the Cp′ ligand could complicate the interactions between the boron cocatalyst and the catalyst.

Structure-properties relationship
To find any factors that could account for the activation barriers of those five catalysts shown in Fig. 4c, we analyzed several electronic and geometric parameters, including HOMO-LUMO energies, atomic charge, and distances (M-Cen and C α -O α ) involved in the catalytic center and the results are collected in Table 3.The Mulliken charges (q) of the metal center were analyzed at the Complex species and we found that the positive charge on the Zr atom correlates well with the observed activity difference of the Cp′ catalysts (C1, C3, and H1), i.e. the higher the charge on the Zr metal (q, C1 < C3 < H1), the higher the activation energies (Ea1, C1 < C3 < H1).Similar trends have been reported previously for styrene polymerization by a rare-earth-metal catalyst 45 .This is related to the donor-acceptor interaction between the ligand and the metal 46 as can be seen from the longer distances between the Zr atom and the centroid of the Cp ring (Zr-Cen = 2.246, 2.274, and 2.289 Å in 2-C1, 2-C3, and 2-H1, respectively, Table 3).There are significant trends between the energy gaps and the number of the methyl group on the Cp ring of the catalyst (ΔE values in Table 3) as well as those between the activation energies (Ea1) and the Zr-Cen distances, atomic charge, and ΔE values (Figure S5).
To evaluate further whether the naked cation model could be a representative of the metallocene/borate system, we compared the two activation Gibbs energies of the rate-limiting step among the three catalytic systems  www.nature.com/scientificreports/(Fig. 7).We observed significant correlations between the activation energies, Ea1, from the naked model and the full models of [B(C 6 F 5 ) 4 ] -and [MeB(C 6 F 5 ) 3 ] -containing systems, with R 2 of 0.6468 and 0.6905, respectively.These relationship analyses further suggest us that it is possible to use the naked cation approach as a pre-screening strategy for catalyst design for CROP of cyclic ester.

Catalyst design with the naked model
Previous structure-activity relationship studies have proven that structural and electronic properties of metallocene catalyst can be tunable by quantum chemistry calculations for CROP [17][18][19]36 and olefin polymerization 27,28,47,48 .
As the results shown above for a small set of zirconocenes, the naked cation approach could be used as a simple strategy in pre-screening a large data set of compounds for candidate catalysts (Fig. 8).To expand the data set, we employed a larger data set of 54 zirconocene structures and used them for in silico pre-screening calculations with the naked model.We also used the same energetic parameters previously described for the CROP activity, i.e., ΔG com for the relative complex stability (Cation → Complex) and Ea1 for the activation energy in the first propagation step (R1 → TS1).From this data set, 29 optimized structures were taken from our previous study 16 and additional 25 structures were built and optimized, leading to a total of 54 catalyst structures.Nomenclature of each catalyst structures (C n and H n ) was assigned in a similar fashion as previously reported by Karttunen et al. 48where C n and H n represent the available crystallographic structures and in silico design, respectively.The resulting structures and energetics as well as electronic properties for the 54 zirconocenes are collected in   Tables S2 and S3.Bar charts showing the ∆G com and Ea1 values are plotted in Fig. 9. Using the C1 energies as a reference, it is clear that both ancillary Cp′ ligand and ligand/fluorine substitutions significantly alter the ∆G com and Ea1 values, and the latter has a strong effect.Other structural and electronic characteristics of zirconocene catalysts in CL polymerizations were described as ancillary ligand, ligand substituent and bridge below.

Effect of ancillary Cp′ ligand
The effect of Cp′ ligands on electronic and steric properties of the complexes can be seen by considering the three Cp′ ligands (Cp, Ind, Flu).Given the data in Tables S2 and S3 and Fig. 9, it is seen that changing Cp into more electron-rich ligand Ind or Flu clearly affects both electronic and steric nature of the catalyst in the process.This type of modification has a greater effect on the complex formation (up to 29 kJ/mol) as it increases the complex stability with a trend C1 ≈ C15 < H16 < H22 (corresponding to the ∆G values of − 98.0, − 97.9, − 122.3 − 127.3 kJ/mol, respectively, Table S2).However, this effect yields a slight increase in the Ea1 values in the order: C1 (74.5 kJ/mol) < C15 (77.4 kJ/mol) < H16 (80.7 kJ/mol) < H22 (80.8 kJ/mol).These results are in-line with our previous study 17 and support the slower kinetics observed for larger ancillary ligands 24 .
The changes in electronic properties due to the larger ancillary ligands can be seen from the atomic charge and HOMO-LUMO energies.The positive charges of Zr atoms become higher as the reaction proceeds from Cation to Complex (Cation/Complex): C1 (1.537e − /1.861e − ), C15 (1.807e − /1.893e − ), H16 (1.672e − /2.007e − ), and H22 (1.647e − /1.886e − ).The HOMO energies are strongly perturbed, with the Cation species being most affected (Figure S6a).This ancillary effect also shows a trend with ∆E HOMO-LUMO values at Complex and R1 (Figure S6b).The electronic origin of this modification is evident in the complexation energies of H16 and H17 (− 122.3 vs − 85.1 kJ/mol, respectively).The loss of a benzene ring in H16 into a cyclohexane in H17 creates a steric clash www.nature.com/scientificreports/ between the cyclohexane ring and the methylene carbon of the CL monomer (Figure S7).However, this effect is insignificant to the activation barrier of the first propagation (80.7 vs 79.9 kJ/mol for H16 and H17).

Effect of bridge
The effect of bridge can be seen by comparing C1 to C5, C6, H4, H6 and C14, C3 to C7 and C8, H1 to C10, C11, C12, and C13, H16 to H18, C22, and C18, H22 to C24 and H23, and others (H2 vs H3, C15 vs C16, H7 vs H8, H10 vs H13, and H24 vs H25).As in other ligand effects, the bridge ligands generally destabilize all intermediates, leading to the lower stabilization to the complex and the higher activation energies to the first propagation.Longer bridge ligands also show similar effects, as can be seen in C1, C5, and C14 for Cp 2 ligand or H16, C22, C18 for Ind 2 ligand (Figure S18).This can be explained from the structural constraint found along the reaction intermediates, which can be justified from the Cen-Zr-Cen (α) and Cen-Cen (β) angles where Cen is the centroid of the Cp′ ring.We analyzed these angles for all 54 catalysts and found three different clusters (Fig. 10): carbon bridge ligands (8 catalysts; cluster I), non-carbon bridge ligand (Si, Ge, P; 21 catalysts; cluster II) and unbridged ligands (25 catalysts; cluster III).All these data can be found in Tables S2 and S3.The cluster I show more open structures of catalyst, which was characterized by the highest β values of around 69°-75°.This can be rationalized from the smallest atomic radius of bridging carbon compared with the other atoms (Si, Ge, P).Cluster II also shows similar characteristics of the bridged ligand but to a lesser extent (β ~ 55° to 65° and α ~ 122° to 132°).The cluster III is the largest members representing the unbridged ligands, with a narrow Cp-Cp structure (β ~ 34° to 55° and α ~ 131° to 142°).Excellent linear correlation (R 2 = 0.91668) is found when plotting between the α and β angles at Cation and the values are significantly reduced when going to Complex and R1 (Fig. 10).
Furthermore, bridge ligands can be used for selective coordination of monomers, particularly for CpInd ligands.They can be used for selective insertion from either a five-membered or a six-membered ring.Comparing C15 (unbridged) with C16 (bridged), the unbridged case shows similar energetics for both insertions C15(5) and C15(6) (− 98.0 vs − 97.9 kJ/mol for ∆G com and 76.1 vs 77.4 kJ/mol for Ea1, respectively).However, due to the bridge, complex formations with the five-membered insertion of C16(5) are 26 kJ/mol more stable than those with the six-membered insertion of C16(6), while the activation energy remains unchanged (76.0 vs 74.1 kJ/mol for C16(5) and C16(6), Figure S19).
Despite the destabilizing effect by the bridge, there is an exception when the Cp′ ligand substituted by the electron-donating group.This can be seen in H3, which possess butyl substituent on the Cp ring and a silicon bridged atom, is found to exhibit the lower activation energy (69.9 kJ/mol) compared to that of C1 (Fig. 9b).Overall, the results indicate that the catalytic performance relies not only on the catalyst structure but also on the electronic nature of the ligand, which has a profound effect on the reaction performance.

Conclusion
We present a structure-property relationship in the CL CROP reaction on the 54 zirconocene structures using DFT calculations.Catalyst structure effect is then investigated as an important factor for tuning the CROP activity, which is analyzed in terms of ancillary ligands, ligand substituents, and bridging units.The polymerization was confirmed to proceed via a cationic mechanism where the ring opening taking place via alkyl-bond cleavage.The activation free energies calculated for the insertion of the first monomer into the C α -O α bond within the CL ring for C1, C3, and H1 does correlate with the observed difference in the catalyst activity reported experimentally.The calculations also show that ligand substituent has a dominant effect to the reaction, especially the electron-donating groups of alkyl and phenyl substituents.The ancillary ligands showed a direct effect on electronic energies, while the bridges are mainly responsible for the catalyst geometry leading to a decrease in the complex stability and an increase in the activation energy of the propagating step.From the structure-property analysis, we demonstrate that the naked model approach can be applied to prescreen catalyst candidate for the CROP process.Overall, this study contributes to a better understanding of the structural characteristics of metallocene catalysts in CL polymerizations, paving the way for further investigations and advancements in this field.

Figure 2 .
Figure 2. Mechanism of CROP of CL using zirconocene/borate systems.Important species involved in the initiation and propagation mechanism, including isolated Cp 2 ZrMe + and CL, monomer-activated complex, reactant, transition state and product, are also shown and labeled as Cation, Complex, R n , TS n and P n , respectively.

Figure 4 .Figure 5 .
Figure 4. Heatmap corresponding to the energy cross-correlation between the catalyst pair for the (a) initiation and (b) first propagation.(c) Bar chart showing activation free energies (ΔG ǂ Ea1 , kJ/mol) of the first propagation step using the selected catalysts (C1, C3, H1, C4 and H16) for the three catalytic systems (C + , 1, and 2).

Figure 6 .
Figure 6.Optimized transition-state structures (2-TS1) in the first propagation step for different [X][B(C 6 F 5 ) 4 ] - (2) systems (X = C1, C3, H1, C4 and H16).Atoms in the top view are shown in the ball-and-stick model.All distances are in units of Å. Imaginary frequencies for each TS are also indicated.Atomic definitions discussed in the text are also given in the figure.

Figure 7 .
Figure 7. Relationship between the two activation Gibbs energies (ΔG ǂ Ea1 ) obtained from different catalytic model systems: (a) naked model vs B(C 6 F 5 ) 4 ; (b) naked model vs MeB(C 6 F 5 ) 3 and (c) MeB(C 6 F 5 ) 3 vs B(C 6 F 5 ) 4.Values are in units of kJ/mol.C1 is an outliner, which does not include in the equation.

Figure 9 .
Figure 9. Bar chart showing (a) complexation energies, ∆G com , and (b) activation free energies of the first propagation (Ea1) for different catalyst ligands.The energies of C1 (∆G com = − 98.0 kJ/mol and Ea1 = 74.5 kJ/ mol) are indicated with dashed lines.

Table 3 .
Calculated Mulliken atomic charge (q, in unit of e) on Zr atom and HOMO and LUMO energies (E HOMO and E LUMO , au) a of Complex and energy gaps (ΔE, au) a between the HOMO of CL and the LUMO of Complex for the naked cation (C + ) and [X][A¯] systems (X = C1